Analiza zbioru danych Bodyfat
Słowem wstępu
Dane pochodzą ze strony internatowej Kaggle.com ze zbioru danych Bodyfat.
Dane wymieniają szacunkowe wartości procentowe tkanki tłuszczowej mierzone pod wodą, wagę i pomiary różnych obwodów ciała 252 mężczyzn.
Rows: 252 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (15): Density, BodyFat, Age, Weight, Height, Neck, Chest, Abdomen, Hip, ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
| Density | BodyFat | Age | Weight | Height | Neck | Chest | Abdomen | Hip | Thigh | Knee | Ankle | Biceps | Forearm | Wrist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.07 | 12.3 | 23 | 154.25 | 67.75 | 36.2 | 93.1 | 85.2 | 94.5 | 59.0 | 37.3 | 21.9 | 32.0 | 27.4 | 17.1 |
| 1.09 | 6.1 | 22 | 173.25 | 72.25 | 38.5 | 93.6 | 83.0 | 98.7 | 58.7 | 37.3 | 23.4 | 30.5 | 28.9 | 18.2 |
| 1.04 | 25.3 | 22 | 154.00 | 66.25 | 34.0 | 95.8 | 87.9 | 99.2 | 59.6 | 38.9 | 24.0 | 28.8 | 25.2 | 16.6 |
| 1.08 | 10.4 | 26 | 184.75 | 72.25 | 37.4 | 101.8 | 86.4 | 101.2 | 60.1 | 37.3 | 22.8 | 32.4 | 29.4 | 18.2 |
| 1.03 | 28.7 | 24 | 184.25 | 71.25 | 34.4 | 97.3 | 100.0 | 101.9 | 63.2 | 42.2 | 24.0 | 32.2 | 27.7 | 17.7 |
| 1.05 | 20.9 | 24 | 210.25 | 74.75 | 39.0 | 104.5 | 94.4 | 107.8 | 66.0 | 42.0 | 25.6 | 35.7 | 30.6 | 18.8 |
Objaśnienie zmiennych zbioru danych:
Density- Gęstość wyznaczona na podstawie ważenia pod wodąBodyFat- Procent tkanki tłuszczowej z równania Siri (1956)Age- Wiek (lata)Weight- Waga (w funtach)Height- Wysokość (cale)Neck- Obwód szyi (cm)Chest- Obwód klatki piersiowej (cm)Abdomen- Obwód brzucha (cm)Hip- Obwód bioder (cm)Thigh- Obwód uda (cm)Knee- Obwód kolana (cm)Ankle- Obwód kostki (cm)Biceps- Obwód bicepsa (cm)Forearm- Obwód przedramienia (cm)Wrist- Obwód nadgarstka (cm)
\hspace{0.8cm} Fat(\%) = \frac{495}{density}-450
Przekształcimy wagę z futów na kilogramy i wysokość w calach na cm według standardu środkowo-europejskigo.
Code
dane["Weight2"] <- round(dane$Weight*0.435,2)
dane["Height2"] <- round(dane$Height*2.54,2)Dodamy kolumnę ID dla łatwiejszej obróbki danych.
Code
ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)Również dodamy nową kolumnę BMI - to wskaźnik masy ciała, który każdy jest w stanie wyznaczyć i zinterpretować w której kategorii znajduje się jego wynik: niedowaga, waga prawidłowa, nadwaga i otyłość. Uważamy, że nowa zmienna BMI będzie dobrym dopełnieniem naszego zbioru danych.
Code
dane["BMI"] <- round(dane$Weight2/(dane$Height2/100)^2)Dodamy kolumnę category dla zmiennej BMI, która ma 3 stany:
- za niska
- w normie
- za wysoka
Code
dane["Category"] <- as.factor(ifelse(dane$BMI< 20.1, '1',
ifelse(dane$BMI< 25.9, '2','3')))Dodamy też kategorie dla zmiennej BodyFat.
Code
dane["Category2"] <- as.factor(ifelse(dane$BodyFat< 11, '1',
ifelse(dane$BodyFat< 24, '2','3')))Sprawdzimy czy mamy braki danych w poszczególnych kolumnach:
Code
apply(dane, 2, function(x) sum(is.na(x))) ID Density BodyFat Age Weight Height Neck Chest
0 0 0 0 0 0 0 0
Abdomen Hip Thigh Knee Ankle Biceps Forearm Wrist
0 0 0 0 0 0 0 0
Weight2 Height2 BMI Category Category2
0 0 0 0 0
Nie mamy braków danych w zbiorze. Kolejnym krokiem będzie wyznaczenie podstawowych statystyk (bez kolumny ID):
Code
summary(dane[,-1]) Density BodyFat Age Weight
Min. :0.995 Min. : 0.00 Min. :22.00 Min. :118.5
1st Qu.:1.041 1st Qu.:12.47 1st Qu.:35.75 1st Qu.:159.0
Median :1.055 Median :19.20 Median :43.00 Median :176.5
Mean :1.056 Mean :19.15 Mean :44.88 Mean :178.9
3rd Qu.:1.070 3rd Qu.:25.30 3rd Qu.:54.00 3rd Qu.:197.0
Max. :1.109 Max. :47.50 Max. :81.00 Max. :363.1
Height Neck Chest Abdomen
Min. :29.50 Min. :31.10 Min. : 79.30 Min. : 69.40
1st Qu.:68.25 1st Qu.:36.40 1st Qu.: 94.35 1st Qu.: 84.58
Median :70.00 Median :38.00 Median : 99.65 Median : 90.95
Mean :70.15 Mean :37.99 Mean :100.82 Mean : 92.56
3rd Qu.:72.25 3rd Qu.:39.42 3rd Qu.:105.38 3rd Qu.: 99.33
Max. :77.75 Max. :51.20 Max. :136.20 Max. :148.10
Hip Thigh Knee Ankle Biceps
Min. : 85.0 Min. :47.20 Min. :33.00 Min. :19.1 Min. :24.80
1st Qu.: 95.5 1st Qu.:56.00 1st Qu.:36.98 1st Qu.:22.0 1st Qu.:30.20
Median : 99.3 Median :59.00 Median :38.50 Median :22.8 Median :32.05
Mean : 99.9 Mean :59.41 Mean :38.59 Mean :23.1 Mean :32.27
3rd Qu.:103.5 3rd Qu.:62.35 3rd Qu.:39.92 3rd Qu.:24.0 3rd Qu.:34.33
Max. :147.7 Max. :87.30 Max. :49.10 Max. :33.9 Max. :45.00
Forearm Wrist Weight2 Height2
Min. :21.00 Min. :15.80 Min. : 51.55 Min. : 74.93
1st Qu.:27.30 1st Qu.:17.60 1st Qu.: 69.16 1st Qu.:173.35
Median :28.70 Median :18.30 Median : 76.78 Median :177.80
Mean :28.66 Mean :18.23 Mean : 77.83 Mean :178.18
3rd Qu.:30.00 3rd Qu.:18.80 3rd Qu.: 85.69 3rd Qu.:183.52
Max. :34.90 Max. :21.40 Max. :157.97 Max. :197.49
BMI Category Category2
Min. : 17.00 1: 29 1: 48
1st Qu.: 22.00 2:145 2:128
Median : 24.00 3: 78 3: 76
Mean : 24.87
3rd Qu.: 26.00
Max. :159.00
Zauważyłyśmy, że wartość minimalna zmiennej BodyFat jest równa 0. Przyjrzymy się tej obserwacji.
Code
dane %>%
filter(BodyFat == 0) ID Density BodyFat Age Weight Height Neck Chest Abdomen Hip Thigh Knee Ankle
1 182 1.1089 0 40 118.5 68 33.8 79.3 69.4 85 47.2 33.5 20.2
Biceps Forearm Wrist Weight2 Height2 BMI Category Category2
1 27.7 24.6 16.5 51.55 172.72 17 1 1
Przypuszczamy, że to może być błąd we wprowadzeniu danych.
Code
495/(dane[182,"Density"])-450[1] -3.611687
Obliczając poziom tłuszczu z równania Siri uzyskałyśmy wartość ujemną, co oznacza, że wzór Siri nie jest optymalny dla każdego przypadku. Ze względu na możliwą do uzyskania wysokość i wagę zostawiamy obserwację w zbiorze danych ze zmienionym poziomem tłuszczu.
Code
dane[182,2] <- 495/(dane[182,"Density"])-450Przedstawimy podstawowe statystyki danych w bardziej zgrabnej formie:
Code
st_op <- round(apply(dane[,c(1,2,3,8,16,17,18)],2,summary),2)
st_op <- rbind(st_op,St.dev=apply(dane[,c(1,2,3,8,16,17)],2,sd))
colnames(st_op) <- c("Density", "BodyFat","Age" ,"Abdomen", "Weight", "Height", "BMI")
rownames(st_op) <- c('minimum','kwantyl dolny','mediana','średnia','kwantyl górny','maksimum','odchylenie standardowe')
as.data.frame(round(st_op,2)) %>%
kable(caption="Statystyki zbioru")| Density | BodyFat | Age | Abdomen | Weight | Height | BMI | |
|---|---|---|---|---|---|---|---|
| minimum | 1.00 | -3.61 | 0.00 | 79.30 | 15.80 | 51.55 | 74.93 |
| kwantyl dolny | 63.75 | 1.04 | 12.47 | 94.35 | 17.60 | 69.16 | 173.35 |
| mediana | 126.50 | 1.05 | 19.20 | 99.65 | 18.30 | 76.78 | 177.80 |
| średnia | 126.50 | 1.04 | 19.15 | 100.82 | 18.23 | 77.83 | 178.18 |
| kwantyl górny | 189.25 | 1.07 | 25.30 | 105.38 | 18.80 | 85.69 | 183.52 |
| maksimum | 252.00 | 1.10 | 47.50 | 136.20 | 21.40 | 157.97 | 197.49 |
| odchylenie standardowe | 72.89 | 0.29 | 8.37 | 8.43 | 0.93 | 12.78 | 72.89 |
Przedstawimy również graficzną prezentację rozkładów :
Code
ggplot(dane, aes(x=Weight2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(50, 160, 10)) +
labs(title = "Histogram wagi [kg]", x = "Waga [kg]", y = "Liczebność") + theme_bw()Code
ggplot(dane, aes(x=Height2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(74, 204, 10)) +
labs(title = "Histogram wzrostu [cm]", x = "Wzrost [cm]", y = "Liczebność") + theme_bw()Na wykresach widzimy, że niektóre obserwacje znacznie wychodzą poza tendencję, w dalszym ciągu zobaczymy jak to wpłynie na wyniki.
Ze statystyk opisowych odczytałyśmy wartość maksymalną BMI, która wyniosła 159. Wyświetlimy obserwację posiadającą niemożliwy poziom wskaźnika masy ciała:
Code
dane[which(dane$BMI == 159), c("ID","Weight2", "Height2", "Age")] ID Weight2 Height2 Age
42 42 89.17 74.93 44
Jest ona najprawdopodobniej pomyłką we wprowadzaniu danych, przejrzymy się tej obserwacji w dalszych analizach.
Na koniec przedstawimy wykresy, który ilustrują podział wskaźnika masy ciała i poziomu tłuszczu na kategorie: za niska, w normie, za wysoka.
Code
dane[,-1] %>%
ggplot(aes(x = Category))+
geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
labs(title = "Rozkład BMI", x = "kategoria wagi", y = "liczebność")+
scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
guides(fill=guide_legend(title="kategoria"))+
scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
theme_minimal()Najbardziej liczną jest grupa z prawidłową masą ciała. Tylko 29 osób z badanych 252 mieści się w kategorii BMI poniżej normy.
Code
dane[,-1] %>%
ggplot(aes(x = Category2))+
geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
labs(title = "Rozkład BodyFat", x = "kategoria poziomu tłuszczu", y = "liczebność")+
scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
guides(fill=guide_legend(title="kategoria"))+
scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
theme_minimal()Rozkład poziomu tłuszczu w grupie badanych mężczyzn zachowuje się podobnie jak rozkład BMI.
Code
dane <- dane[,-21]Hipotezy
Czy poszczególne obwody: szyja, klatka piersiowa, brzuch, biodra, uda, kolana, kostka, biceps, przedramię i nadgarstek mają wpływ na BMI?
Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?
Funkcje
Code
gauss_markov <- function(mod) {
tabelka <- NULL
wiersz <- NULL
bp <- bptest(mod)
gq <- gqtest(mod)
dw <- dwtest(mod)
bg <- bgtest(mod, order.by = ~fitted(mod), order = 3)
sh <- shapiro.test(mod$residuals)
ks <- ks.test(mod$residuals, "pnorm")
lil <- nortest::lillie.test(mod$residuals)
res <- resettest(mod, type = "regressor")
rai <- raintest(mod)
har <- harvtest(mod)
w <- c(bp$p.value, gq$p.value, dw$p.value, bg$p.value, sh$p.value, ks$p.value, lil$p.value, res$p.value, rai$p.value, har$p.value)
h0 <- c("jednorodność wariancji", "jednorodność wariancji", "brak autokorelacji reszt","brak autokorelacji do rzędu 3","rozklad reszt jest normalny", "rozklad reszt jest normalny", "rozklad reszt jest normalny", "zależność jest liniowa", "zależność jest liniowa", "zależność jest liniowa")
h1 <- c("brak jednorodności wariancji", "brak jednorodności wariancji", "występuje autokorelacja reszt", "występuje autokorelacja rzędu 3", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "brak zależności liniowej", "brak zależności liniowej", "brak zależności liniowej")
n <- c("Test Breuscha-Pagana", "Test Goldfelda-Quandta", "Test Durbina-Watsona", "Test Breuscha-Godfreya", "Test Shapiro-Wilka", "Test Kolmogorova-Smirnova", "Test Lillieforsa", "Test RESET Ramseya", "Test Rainbow Uttsa", "Test Harveya-Colliera")
for (i in 1:length(w)) {
wiersz <- c(n[i], round(w[i], 3),
if (w[i]<0.05) {
h1[i]} else {h0[i]} )
tabelka <- rbind(tabelka, wiersz)
}
colnames(tabelka) <- c("Nazwa testu", "P-wartość", "Wniosek")
rownames(tabelka) <- c(1:length(n))
tabelka %>%
kable(caption="Założenia Gaussa Markowa")
}Code
odstajace <- function(mod, dane, stopien) {
h <- hatvalues(mod)[hatvalues(mod) > 2*3/nrow(dane)] %>% as.data.frame()
h <- rownames(h)
rsr <- rstandard(mod)[abs(rstandard(mod)) > 2] %>% as.data.frame()
rsr <- rownames(rsr)
rst <- rstudent(mod)[abs(rstudent(mod)) > 3] %>% as.data.frame()
rst <- rownames(rst)
c <- cooks.distance(mod)[cooks.distance(mod) > 4/nrow(dane)] %>% as.data.frame()
c <- rownames(c)
w <- c(h, rsr, rst, c)
b <- unique(w)
counts <- NULL
for (i in 1:length(b)) {
counts <- c(counts, sum(w==b[i]))
}
counts <- cbind(b, counts)
a <- c("obserwacja", "ile spełnia kryteriów")
counts <- counts[which(counts[,2]>=stopien),]
colnames(counts) <- a
counts %>% kable(caption="Obserwacje odstające")
}Code
odstajace_col <- function(dat){
a <- c()
for(i in c(7:16, 19)) {
a <- c(a, dat %>% identify_outliers(colnames(dat)[i]) %>% select(ID))
}
names(a) <- colnames(dat)[c(7:16, 19)]
print(a)
a <- plyr::ldply(a, data.frame)
repetitions <- names(which(table(a[,2])>=5))
print("Obserwacje występujące co najmniej 5 razy: ")
print(repetitions)
}Hipoteza 1
BMI to inaczej wskaźnik masy ciała - z angielskiego body mass index - jest to współczynnik powstały przez podzielenie masy ciała podanej w kilogramach przez kwadrat wysokości podanej w metrach. Klasyfikacja wskaźnika BMI została opracowana wyłącznie dla dorosłych.
BMI = \frac{masa[kg]}{wysokość^2[m]}
Code
knitr::include_graphics("bmi.jpg")Model pierwszy
\hat{BMI} = \beta_{0} + \beta_{1} \cdot Neck + \beta_{2} \cdot Chest + \beta_{3} \cdot Abdomen + \beta_{4} \cdot Hip + \beta_{5} \cdot Thigh + \beta_{6} \cdot Knee + \beta_{7} \cdot Ankle + \beta_{8} \cdot Biceps + \beta_{9} \cdot Forearm + \beta{10} \cdot Wrist
Code
model <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane)
summary(model)
Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee +
Ankle + Biceps + Forearm + Wrist, data = dane)
Residuals:
Min 1Q Median 3Q Max
-6.252 -2.111 -0.352 1.163 120.578
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -11.04020 12.45915 -0.886 0.3764
Neck -0.56573 0.41818 -1.353 0.1774
Chest 0.11726 0.17128 0.685 0.4943
Abdomen 0.06517 0.14352 0.454 0.6501
Hip 0.48499 0.22593 2.147 0.0328 *
Thigh 0.26534 0.25046 1.059 0.2905
Knee 0.03295 0.42583 0.077 0.9384
Ankle -0.13217 0.40735 -0.324 0.7459
Biceps -0.10045 0.31902 -0.315 0.7531
Forearm 0.03633 0.37193 0.098 0.9223
Wrist -1.13466 0.94514 -1.201 0.2311
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.163 on 241 degrees of freedom
Multiple R-squared: 0.241, Adjusted R-squared: 0.2095
F-statistic: 7.653 on 10 and 241 DF, p-value: 1.323e-10
Code
summary(dane$BMI-model$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-6.252 -2.111 -0.352 0.000 1.163 120.578
Mediana jest mniejsza od średniej to rozkład reszt jest prawostronnie asymetryczny.
Wykresy diagnostyczne:
Code
autoplot(model)Założenia Gaussa Markowa
Code
gauss_markov(model)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.039 | brak jednorodności wariancji |
| Test Goldfelda-Quandta | 1 | jednorodność wariancji |
| Test Durbina-Watsona | 0.769 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.794 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0 | rozklad reszt nie jest normalny |
| Test Kolmogorova-Smirnova | 0 | rozklad reszt nie jest normalny |
| Test Lillieforsa | 0 | rozklad reszt nie jest normalny |
| Test RESET Ramseya | 0.034 | brak zależności liniowej |
| Test Rainbow Uttsa | 0 | brak zależności liniowej |
| Test Harveya-Colliera | 0.128 | zależność jest liniowa |
Budując model pełny otrzymałyśmy skorygowany R^2 równy 0.21, przyjrzałyśmy się wykresom i zauważyłyśmy, że obserwacja 42 jest nietypowa. Sprawdzimy, czy tak jest na prawdę.
Code
odstajace(model, dane = dane, 2)| obserwacja | ile spełnia kryteriów |
|---|---|
| 42 | 4 |
| 86 | 2 |
Z powyższego zestawienia widać, że obserwacja 42 może wpływać na to, że nasz model nie spełnia założeń. Sprawdźmy czy w wyniku usuwania jej ze zbioru, nasz model się poprawi.
Model drugi
Code
dane2 <- dane[-42,]
model2 <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane2)
summary(model2)
Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee +
Ankle + Biceps + Forearm + Wrist, data = dane2)
Residuals:
Min 1Q Median 3Q Max
-3.0303 -0.6120 -0.0121 0.6024 4.0308
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -12.78573 1.66370 -7.685 3.90e-13 ***
Neck 0.07917 0.05612 1.411 0.159594
Chest 0.13840 0.02287 6.051 5.49e-09 ***
Abdomen 0.11926 0.01917 6.222 2.17e-09 ***
Hip 0.08192 0.03037 2.697 0.007485 **
Thigh 0.12382 0.03347 3.700 0.000268 ***
Knee -0.24602 0.05691 -4.323 2.26e-05 ***
Ankle 0.07697 0.05442 1.414 0.158566
Biceps 0.06284 0.04262 1.474 0.141674
Forearm -0.02033 0.04967 -0.409 0.682587
Wrist -0.00730 0.12658 -0.058 0.954061
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.09 on 240 degrees of freedom
Multiple R-squared: 0.9078, Adjusted R-squared: 0.904
F-statistic: 236.4 on 10 and 240 DF, p-value: < 2.2e-16
Code
summary(dane2$BMI-model2$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.03026 -0.61203 -0.01207 0.00000 0.60240 4.03083
Mediana jest bardzo bliska 0, co wskazuję, na brak asymetrii.
Sprawdzimy jak wyglądają wykresy diagnostyczne:
Code
autoplot(model2)Z wykresów diagnostycznych raczej spodziewamy się niejednorodności wariancji i braku normalności rozkładu reszt. Sprawdzamy to za pomocą testów.
Założenia Gaussa Markowa
Code
gauss_markov(model2)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.001 | brak jednorodności wariancji |
| Test Goldfelda-Quandta | 0.753 | jednorodność wariancji |
| Test Durbina-Watsona | 0.135 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.575 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.024 | rozklad reszt nie jest normalny |
| Test Kolmogorova-Smirnova | 0.814 | rozklad reszt jest normalny |
| Test Lillieforsa | 0.15 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0 | brak zależności liniowej |
| Test Rainbow Uttsa | 0.094 | zależność jest liniowa |
| Test Harveya-Colliera | 0.811 | zależność jest liniowa |
Po sprawdzeniu założeń dochodzimy do wniosku, że jednak większa część z nich na razie nie jest spełniona. Sprawdzimy obserwacje odstające dla poszczególnych zmiennych:
Code
odstajace_col(dane)$Neck
[1] 39 45 106
$Chest
[1] 39 41
$Abdomen
[1] 39 41 216
$Hip
[1] 35 39 41
$Thigh
[1] 39 41 152 169
$Knee
[1] 39 192 244
$Ankle
[1] 31 39 86
$Biceps
[1] 39
$Forearm
[1] 45 159 175 206 226
$Wrist
[1] 39 41 226 252
$BMI
[1] 39 41 42 216
[1] "Obserwacje występujące co najmniej 5 razy: "
[1] "39" "41"
Zauważyłyśmy, że obserwacje 39, 41 są odstające prawie dla każdej zmiennej. Również obserwacja 39 jest w każdej statystyce rozpatrywanej w funkcji odstajace(). Sprawdzimy jak model będzie wygłądać bez tych obserwacji:
Code
dane3 <- dane[-c(39,41,42),]Po wyczyszczeniu danych rozdzielamy nasz zbior na testowy i uczący.
Code
ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)
set.seed(221)
wektor <- sample(1:249, 180, replace = FALSE)
uczacy <- dane3[wektor,]
testowy <- dane3[-wektor,]Dalsze analizy będziemy przeprowadzać na zbiorze uczącym.
Model trzeci
Code
model3 <- lm(BMI~Chest+Neck+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm, data=uczacy)
summary(model3)
Call:
lm(formula = BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee +
Ankle + Biceps + Forearm, data = uczacy)
Residuals:
Min 1Q Median 3Q Max
-2.6920 -0.6510 0.0483 0.6423 3.4233
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.163994 1.854714 -4.941 1.85e-06 ***
Chest 0.139161 0.025767 5.401 2.21e-07 ***
Neck -0.014175 0.059800 -0.237 0.8129
Abdomen 0.137760 0.022117 6.229 3.57e-09 ***
Hip -0.002438 0.036103 -0.068 0.9462
Thigh 0.145343 0.035583 4.085 6.79e-05 ***
Knee -0.169290 0.071710 -2.361 0.0194 *
Ankle 0.020266 0.066001 0.307 0.7592
Biceps 0.042724 0.046566 0.917 0.3602
Forearm 0.122550 0.061477 1.993 0.0478 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.043 on 170 degrees of freedom
Multiple R-squared: 0.8924, Adjusted R-squared: 0.8867
F-statistic: 156.6 on 9 and 170 DF, p-value: < 2.2e-16
Możemy zauważyć, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę.
Code
summary(dane3$BMI-model3$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-11.81280 -3.07204 0.14458 0.06521 2.98281 10.48664
Code
autoplot(model3)Wykresy diagnostyczne wyglądają lepiej, niż w poprzednim modelu, spodziewamy się, że większość zalożeń będzie spełnionych.
Założenia Gaussa Markowa
Code
gauss_markov(model3)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.682 | jednorodność wariancji |
| Test Goldfelda-Quandta | 0.933 | jednorodność wariancji |
| Test Durbina-Watsona | 0.787 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.197 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.867 | rozklad reszt jest normalny |
| Test Kolmogorova-Smirnova | 0.98 | rozklad reszt jest normalny |
| Test Lillieforsa | 0.843 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0.157 | zależność jest liniowa |
| Test Rainbow Uttsa | 0.858 | zależność jest liniowa |
| Test Harveya-Colliera | 0.908 | zależność jest liniowa |
W kolejnym kroku po zbudowaniu pełnego modelu przejrzymy się czy da się nasz model ulepszyć, widzimy, że model posiada 9 zmiennych, co jest sporo, chcemy zredukować je ilość za pomocą metody PCA, którą przedstawiamy poniżej:
Najpierw sprawdzamy korelację zmiennych:
Widać, że niektóre zmienne są mocno skorelowane, więc PCA będzie tutaj dobrym wyborem:
Code
model3.pca <- prcomp(uczacy[,-c(1:6,16:20)], scale = T)
summary(model3.pca)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.5164 0.86709 0.7794 0.6490 0.53794 0.49801 0.43170
Proportion of Variance 0.7036 0.08354 0.0675 0.0468 0.03215 0.02756 0.02071
Cumulative Proportion 0.7036 0.78715 0.8547 0.9014 0.93360 0.96116 0.98187
PC8 PC9
Standard deviation 0.30138 0.26899
Proportion of Variance 0.01009 0.00804
Cumulative Proportion 0.99196 1.00000
Code
fviz_eig(model3.pca, addlabels = TRUE)Code
fviz_contrib(model3.pca, choice = "var", axes = 9)Widzimy, że pierwsze 3 składowych już dałyby nam niemal 86% wyjaśnionej wariancji, co byłoby zgodnie z kryterium.
Tworzymy zmienne PC1,PC2,PC3.
Code
PC1 <- model3.pca$x[,1]
PC2 <- model3.pca$x[,2]
PC3 <- model3.pca$x[,3]Code
model3_2 <- lm(BMI~PC1+PC2+PC3, data=uczacy)
summary(model3_2)
Call:
lm(formula = BMI ~ PC1 + PC2 + PC3, data = uczacy)
Residuals:
Min 1Q Median 3Q Max
-2.8746 -0.7642 0.0386 0.8813 4.0660
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.13333 0.08748 275.868 < 2e-16 ***
PC1 -1.09812 0.03486 -31.500 < 2e-16 ***
PC2 -0.80886 0.10117 -7.995 1.65e-13 ***
PC3 -0.42268 0.11255 -3.756 0.000235 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.174 on 176 degrees of freedom
Multiple R-squared: 0.8588, Adjusted R-squared: 0.8564
F-statistic: 356.8 on 3 and 176 DF, p-value: < 2.2e-16
Code
gauss_markov(model3_2)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.355 | jednorodność wariancji |
| Test Goldfelda-Quandta | 0.949 | jednorodność wariancji |
| Test Durbina-Watsona | 0.872 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.371 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.334 | rozklad reszt jest normalny |
| Test Kolmogorova-Smirnova | 0.362 | rozklad reszt jest normalny |
| Test Lillieforsa | 0.433 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0.03 | brak zależności liniowej |
| Test Rainbow Uttsa | 0.396 | zależność jest liniowa |
| Test Harveya-Colliera | 0.587 | zależność jest liniowa |
R^2 tego modelu to 0.85 Założenia Gaussa-Markowa z różnych testów troszkę się polepszyły.
Model zagnieżdżony
Code
mod_0 <- lm(BMI~1, data = uczacy)
met_tyl <- stats::step(model3, scope = c(mod_0, model3),
direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients(Intercept) Chest Abdomen Thigh Knee Forearm
-9.4910775 0.1444103 0.1339823 0.1549096 -0.1639500 0.1378829
Code
model_zag <- lm(BMI~Chest+Abdomen+Thigh+Knee+Forearm, data=uczacy)
summary(model_zag)
Call:
lm(formula = BMI ~ Chest + Abdomen + Thigh + Knee + Forearm,
data = uczacy)
Residuals:
Min 1Q Median 3Q Max
-2.6718 -0.6431 0.0511 0.6350 3.5310
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.49108 1.51928 -6.247 3.11e-09 ***
Chest 0.14441 0.02415 5.981 1.23e-08 ***
Abdomen 0.13398 0.01978 6.774 1.85e-10 ***
Thigh 0.15491 0.02801 5.531 1.15e-07 ***
Knee -0.16395 0.06147 -2.667 0.00837 **
Forearm 0.13788 0.05417 2.545 0.01178 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.033 on 174 degrees of freedom
Multiple R-squared: 0.8918, Adjusted R-squared: 0.8887
F-statistic: 286.7 on 5 and 174 DF, p-value: < 2.2e-16
Z podsumowania można wnioskować, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę, tak samo jak w modelu pełnym.
Code
summary(uczacy$BMI-model_zag$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.67183 -0.64313 0.05105 0.00000 0.63502 3.53102
Założenia Gaussa Markowa
Code
gauss_markov(model_zag)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.638 | jednorodność wariancji |
| Test Goldfelda-Quandta | 0.918 | jednorodność wariancji |
| Test Durbina-Watsona | 0.787 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.63 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.741 | rozklad reszt jest normalny |
| Test Kolmogorova-Smirnova | 0.979 | rozklad reszt jest normalny |
| Test Lillieforsa | 0.808 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0.277 | zależność jest liniowa |
| Test Rainbow Uttsa | 0.813 | zależność jest liniowa |
| Test Harveya-Colliera | 0.945 | zależność jest liniowa |
Code
autoplot(model_zag)Code
anova(model3, model_zag)Analysis of Variance Table
Model 1: BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee + Ankle + Biceps +
Forearm
Model 2: BMI ~ Chest + Abdomen + Thigh + Knee + Forearm
Res.Df RSS Df Sum of Sq F Pr(>F)
1 170 184.79
2 174 185.82 -4 -1.0382 0.2388 0.9161
Na podstawie analizy wariancji przyjmujemy, że model prostszy jest lepszy.
Predykcja
Code
pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,19)]Code
cbind(pred, fakt = fakt$BMI) %>%
ggplot(aes(x = pred, y = fakt))+
geom_point()+
geom_smooth()`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Podsumowanie
Postać modelu:
\hat {BMI} = -9.49 + 0.14*Chest + 0.13*Abdomen + 0.15*Thigh - 0.16*Knee + 0.14*Forearm
Nasz model wyjaśnia niemal 89% wariancji. Zmienne wchodzące w skład modelu to obwody: klatka piersiowa, brzuch, uda, kolana i przedramię. Wskazuje to na dużą zależność tych miar ze wskaźnikiem masy ciała. Oznacza to, że dla przeciętnych mężczyzn wskaźnik ten jest wystarczająco dobrą miarą aby określić swój stan fizyczny. Pamiętajmy jednak, że przy dużej objętości tkanki mięśniowej nie jest on miarodajny i raczej wybralibyśmy procentowy wskaźnik poziomu tłuszczu.
Hipoteza 2
Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?
Code
ggplot(dane3,aes(x=BodyFat,y=Weight2)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i wagi") +
geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 250)`geom_smooth()` using formula = 'y ~ x'
Code
ggplot(dane3,aes(x=BodyFat,y=Abdomen)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i obwodu brzucha") +
geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 115)`geom_smooth()` using formula = 'y ~ x'
Model pierwszy
\hat{BodyFat} = \beta_{0} + \beta_{1} \cdot Chest + \beta_{2} \cdot Abdomen + \beta_{3} \cdot Hip + \beta_{4} \cdot Thigh + \beta_{5} \cdot Knee + \beta_{6} \cdot Ankle + \beta_{7} \cdot Biceps + \beta_{8} \cdot Forearm + \beta_{9} \cdot Age + \beta{10} \cdot Weight2 + \beta{11} \cdot Height2+\beta{12} \cdot Wrist
Code
model <- lm(BodyFat~Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Age+Weight2+Height2+Wrist, data=uczacy)
summary(model)
Call:
lm(formula = BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee +
Ankle + Biceps + Forearm + Age + Weight2 + Height2 + Wrist,
data = uczacy)
Residuals:
Min 1Q Median 3Q Max
-9.1699 -3.0621 -0.0413 2.8888 9.0910
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -24.67529 25.91811 -0.952 0.34245
Chest -0.05061 0.11979 -0.422 0.67321
Abdomen 0.81385 0.10564 7.704 1.11e-12 ***
Hip -0.10127 0.15618 -0.648 0.51761
Thigh 0.32767 0.16133 2.031 0.04384 *
Knee 0.32390 0.31165 1.039 0.30017
Ankle -0.10546 0.28083 -0.376 0.70774
Biceps 0.21557 0.18828 1.145 0.25388
Forearm 0.37294 0.25023 1.490 0.13801
Age 0.06188 0.03637 1.701 0.09073 .
Weight2 -0.21411 0.17333 -1.235 0.21847
Height2 -0.07795 0.08480 -0.919 0.35931
Wrist -1.95932 0.61351 -3.194 0.00168 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.097 on 167 degrees of freedom
Multiple R-squared: 0.7645, Adjusted R-squared: 0.7476
F-statistic: 45.17 on 12 and 167 DF, p-value: < 2.2e-16
Code
summary(uczacy$BodyFat-model$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.16994 -3.06208 -0.04125 0.00000 2.88881 9.09097
Założenia Gaussa Markowa
Code
gauss_markov(model)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.866 | jednorodność wariancji |
| Test Goldfelda-Quandta | 0.578 | jednorodność wariancji |
| Test Durbina-Watsona | 0.987 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.859 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.226 | rozklad reszt jest normalny |
| Test Kolmogorova-Smirnova | 0 | rozklad reszt nie jest normalny |
| Test Lillieforsa | 0.359 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0.352 | zależność jest liniowa |
| Test Rainbow Uttsa | 0.845 | zależność jest liniowa |
| Test Harveya-Colliera | 0.641 | zależność jest liniowa |
Model zagniezdżony
Code
mod_0 <- lm(BodyFat~1, data = uczacy)
met_tyl <- stats::step(model, scope = c(mod_0, model),
direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients (Intercept) Abdomen Thigh Forearm Age Weight2
-44.07772256 0.83167546 0.41690583 0.51621459 0.07463536 -0.30133107
Wrist
-1.82469765
Code
model_zag <- lm(BodyFat~Abdomen+Thigh+Weight2+Wrist, data=uczacy)
summary(model_zag)
Call:
lm(formula = BodyFat ~ Abdomen + Thigh + Weight2 + Wrist, data = uczacy)
Residuals:
Min 1Q Median 3Q Max
-9.5177 -2.8102 -0.0278 3.1434 9.9275
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -41.85380 9.78803 -4.276 3.12e-05 ***
Abdomen 0.91504 0.06283 14.563 < 2e-16 ***
Thigh 0.34524 0.12056 2.864 0.004701 **
Weight2 -0.31348 0.08241 -3.804 0.000196 ***
Wrist -1.08634 0.52264 -2.079 0.039118 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.143 on 175 degrees of freedom
Multiple R-squared: 0.7476, Adjusted R-squared: 0.7419
F-statistic: 129.6 on 4 and 175 DF, p-value: < 2.2e-16
Możemy zauważyć, że prawdziwy poziom tłuszczu odchyla się od linii regresji średnio o około 4 punkty procentowe.
Code
summary(uczacy$BodyFat-model_zag$fitted.values) Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.51767 -2.81022 -0.02781 0.00000 3.14340 9.92750
Code
autoplot(model_zag)Założenia Gaussa Markowa
Code
gauss_markov(model_zag)| Nazwa testu | P-wartość | Wniosek |
|---|---|---|
| Test Breuscha-Pagana | 0.811 | jednorodność wariancji |
| Test Goldfelda-Quandta | 0.455 | jednorodność wariancji |
| Test Durbina-Watsona | 0.974 | brak autokorelacji reszt |
| Test Breuscha-Godfreya | 0.168 | brak autokorelacji do rzędu 3 |
| Test Shapiro-Wilka | 0.334 | rozklad reszt jest normalny |
| Test Kolmogorova-Smirnova | 0 | rozklad reszt nie jest normalny |
| Test Lillieforsa | 0.25 | rozklad reszt jest normalny |
| Test RESET Ramseya | 0.052 | zależność jest liniowa |
| Test Rainbow Uttsa | 0.885 | zależność jest liniowa |
| Test Harveya-Colliera | 0.718 | zależność jest liniowa |
Zauważamy, że R dopasowane jest zbliżone w modelu zagnieżdżonym do R dopasowanego w modelu pełnym. Oba modele spełniają wszystkie założenia, jedynie test Kołmogorowa - Smirnova wykazuje brak normalności rozkładu reszt. Wszystkie zmienne wchodzące w skład modelu zagnieżdżonego są istotne statystycznie, więc skłanialibyśmy się ku modelowi zagnieżdżonemu.
Code
anova(model, model_zag)Analysis of Variance Table
Model 1: BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee + Ankle + Biceps +
Forearm + Age + Weight2 + Height2 + Wrist
Model 2: BodyFat ~ Abdomen + Thigh + Weight2 + Wrist
Res.Df RSS Df Sum of Sq F Pr(>F)
1 167 2803.8
2 175 3004.4 -8 -200.65 1.4939 0.1628
Na podstawie analizy wariancji stwierdzamy, że model prostszy jest lepszy.
Predykcja
Code
pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,3)]Code
cbind(pred, fakt = fakt$BodyFat) %>%
ggplot(aes(x = pred, y = fakt))+
geom_point()+
geom_smooth()`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Podsumowanie
Postać modelu:
\hat {BodyFat} = -41.85 + 0.92*Abdomen + 0.35*Thigh - 0.31*Weight2 - 1.09*Wrist
Nasz model wyjśnia ponad 74% wariancji. Miary wchodzące w skład modelu to obwody: brzuch, uda, nadgarstek oraz waga w kilogramach, dzięki czemu każdy (mężczyzna) może łatwo sobie wyznaczyć procentowy poziom tłuszczu i sprawdzić do jakiej kategori należy. Pomaga to w kontrolowaniu kondycji fizycznej, która ma wysoki wpływ na zdrowie człowieka.